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Abstract 

We present a fast and efficient tight-binding (TB) method for simulating scanning tunneling 
microscopy (STM) imaging of adsorbate molecules on ultrathin insulating films. Due to the elec- 
tronic decoupling of the molecule from the metal surface caused by the presence of the insulating 
overlayer, the STM images of the frontier molecular orbitals can be simulated using a very effi- 
cient scheme, which also enables the analysis of phase shifts in the STM current. Au-pentacene 
complex adsorbed on a NaCl bilayer on Cu substrate provides an intricate model system, which 
has been previously studied both experimentally and theoretically. Our calculations indicate that 
the complicated shape of the molecular orbitals may cause multivalued constant current surfaces 
— leading to ambiguity of the STM image. The results obtained using the TB method are found 
to be consistent with both DFT calculations and experimental data. 



PACS numbers: 



I. INTRODUCTION 



The continuous miniaturization of electronics to obtain more efficient and smaller com- 
ponents is approaching the classical border beyond which components consisting of single 
molecules are one possible solution. In the development of such molecular components, 
an excellent platform is provided by ultrathin insulating overlayers on metal surfaces com- 
bined with scanning tunneling microscopy (STM). Adsorbates on such films are electroni- 
cally decoupled from the metal substrate, which can be used, e.g., in controlling the charge 
state of the adsorbat^ or studying atoms and molecules without the influence of the metal 
backgrouncP'^. 

For complicated systems theoretical calculations play a critical role in the interpretation 
of the experiments. However, many of the currently studied organic molecules suitable for 
molecular electronic devices are of a scale beyond a feasible system size for ab-initio calcula- 
tions including the substrate. Semiempirical tight-binding (TB) based methods provide an 
alternative solution to this problem, allowing for simulation of larger systems with reason- 
able computational resources. Additionally, TB methods provide different kind of tools to 
analyze both the adsorption and the STM images. 

In this work we study a system consisting of a pentacene molecule bonded with a 
gold atom on a NaCl overlayer previously studied both experimentally and with ab-initio 
methods.- Due to its complicate electronic and geometric structure, this system provides an 
excellent test ground for developing TB based simulations of STM on molecules on insulat- 
ing overlayers. The system exhibits clearly recognizable properties caused by the interaction 
of the pentacene and the gold atom, and by the adsorption of the complex on the insulating 
overlayer, making it ideal for this work. 

We use a purpose-built method for simulating STM in a TB basi^, which has previ- 
ously been applied to studying superconducting materials and organic molecules on metal 
surfaced. This method allows simulations with functionalized tips, analysis of tunneling 
channels, and studying of the molecular orbitals of adsorbed molecules. Here we develop 
the method further by introducing an even faster method for simulating STM of systems 
where a certain molecular orbital dominates the STM image. Provided that there is only 
a single molecular energy level within the bias voltage, the STM current as a function of 
position is approximately proportional to the hopping integral between the STM tip and the 
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frontier molecular orbital. This also allows studying of extremely large systems, or rapidly 
exploring different molecular configurations. The method developed in this work also makes 
it possible to study the phase information of the STM current, which could be utilized, e.g. 
in the simulation of isospectral molecule^. 

This paper is organized as follows: Sec. [ITl explains the theory behind the used calculation 



methods. Sec. [ITT] details the model used for the Au-pentacene calculations. Sec. [TV| presents 
the obtained results, and Sec. |V] discusses the results and their implications. In Appendix 
[A] several different forms of calculating tunneling current are discussed. 

II. THEORY 

We use two different approaches, one based on direct diagonalization of the Hamilto- 
nian and the other on Green's function formalism, to calculate the electronic structure in a 
tight-binding basis. The former approach is used to visualize the molecular orbitals and sim- 
ulate their behaviour in bond formation and adsorption. The latter method allows detailed 
analysis of the electron transport properties and is used in simulating STM images. 

In both cases the Hamiltonian is presented in atomic orbital basis. The on-site terms 
(energies) are parameters fitted to reproduce the density of states compared to density 
functional theory (DFT) calculations, or taken from literature. The off-diagonal elements 



are calculated using modified Slater-Koster hopping integrals described in Sec. [ITT 

In the first method the Hamiltonian is diagonalized using the secular equations He' = 
Ei& in the Hiickel approximation where overlap between different atomic orbitals is ne- 
glected. The obtained eigenvectors c* include weight coefficients = {a\i) of individual 
atomic orbitals a to an eigenstate i of the system with eigenenergy Ei. The weight coeffi- 
cients allow real space projections of the molecular orbitals by using the spatial dependence 
of Slater type orbitals. This makes it also possible to simulate STM with the Tersoff-Hamann 
approaclP (TH) in which surfaces of the constant local density of the sample states at the 
position of the STM tip are identified with the topographic STM images. 

In addition, the contribution of a molecular orbital /i of a free molecule to an electronic 
state 2 of a more complicated system can be obtained by solving the secular equations 
separately for the isolated molecule (non-interacting system) and for the interacting system 
consisting of the molecule and other atoms. Using the eigenvectors and c* calculated in 
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the same basis set a we obtain the weight of molecular orbital [i on the eigenstate i: 

a a 

is actually the eigenvector of Hamiltonian represented in a mixed basis set consisting of 
molecular orbitals of the free molecule and the atomic orbitals of the other atoms of the 
system not involved in the free molecule. 

This formalism can be used both to study the changes in molecular states in the formation 
of chemical bonds, as well as to locate the molecular states of a free molecule from an 
adsorbed system. This also gives information on the interaction of the molecular orbitals 
with the surface states. Since diagonalization of the Hamiltonian is very fast, this method 
allows studying of very large systems. 

In simulating STM images we use a Green's functions approach to calculate the electronic 
structure of the sample, described in detail in Ref. IH The advantage of this approach is ob- 
taining the off-diagonal terms necessary to calculate currents within Todorov-Pendry scheme 
(TP)ISEI. This enables detailed analysis of the STM image; the current can be decomposed to 
the individual contributions of particularly important atomic orbitals which form tunneling 
channel^. As shown in Eq,A4 in Appendix |A[ the current flowing from the STM tip to the 



substrate can be written as 



2er] '"^ 



I{V,n) = Y / jiE,rt)dE (2.2) 







j{E,n) = TT[pl,{E)VtArt)G;,f{E)G%{E)VAn)] 

where denotes the tip position and V the bias voltage. Here j{E,rt) is the differential 
conductance in units of ^ corresponding to bias E/e. Subindices /, t and s refer to final, 
tip and sample states, respectively. and G^ are the retarded and advanced Green's 
function matrices, respectively, and 2// is the full width at half maximum of the broadened 
delta function used in calculating the Green's functions. Vts(rt) = {t{rt)\H\s) is the hopping 
integral from tip to sample state s. The density matrix of the tip p^n{E) is calculated 
separately from the sample. 

On ultrathin insulating films, the STM experiments on adsorbates can produce images 
whose shape is dominated by a single (frontier) orbital of the adsorbate molecul^. In the 
following we show that for this kind of system, STM images can be simulated by hopping 
maps in which the hopping integral from tip to the dominating state of the system is plotted 



as a function of tip position. The idea lies in the decouphng of electronic states of molecule 
and the underlying metal substrate, due to which there is a single eigenstate of the system 
which resembles one of the original molecular orbitals of the adsorbate. Thus, there is a 
subset of sample states (refered to as /i') which correspond to the orbitals of a free molecule 
(refered to as fi). 

First, the tip is approximated by a single electron state so that its density matrix element 



can be taken outside the calculation of matrix trace in Eq. 2.2 In addition, the trace 
can be written double sum 

j{E,n) = p',,{E)J2Y.\^ts{n)G;^{E)\\ (2.3) 

^ / 

Next we diagonalize the Hamiltonian of the sample using the secular equation and rewrite 
the summation keeping in mind that the Vts and are chosen to be real: 

mr,) = pUE)J2J2\^u{E)\'\Gr^{E)\'. (2.4) 

i f 

If in the energy range [E-p, E-p + eV] defined by the bias voltage V there is only a single 
eigenstate z/' of the subset p' of the sample states then we can assume that for all the other 
states 

\VuGif{E)\ < \Vtu'G,,f{E)l for a[\ i ^ v' ^ E e [Ep,Ef + eV] . (2.5) 

The validity of this assumption lies on the nature of the sample states: either the Green's 
function is very small since the state is too far away from the Fermi level (valid for the subset 
fi' and the states corresponding to the atoms of the insulating film), or the exponentially 
decaying hopping integral is negligible due to long inter-atomic distance (valid for the metal 



states). Now the tunneling current of Eq. 2.2 can be simplified to 

iiv,n) ^ '^-^Iv.AnW r pI{e)Y^\g-,^{e)\\ (2.6) 

Jo J 

where the only dependence on the position of the STM tip, rf, is outside of the integral which 
is constant for a certain bias and molecular orbital. Thus, for a bias voltage corresponding 
to a suitable energy range, the spatial dependence of the STM current can be written as: 

i{n) oc = mn)\HW)\' = Y.\^UnY:\\ (2.7) 

a 

In the last form the state v' is represented in atomic orbital basis a. 



Thus, topographic STM images can be approximated with the isosurfaces of the square of 
the hopping integral. However, in the following we shall present results for the bare hopping 
integral as function of tip position — a three-dimensional hopping map. The hopping maps 
reveal the change in the phase of the current, while the constant absolute value surfaces can 
be still compared to topographic STM images. 

To calculate the hopping maps in practise, we only need the hopping integrals between 
the tip and atomic orbitals of the sample, and the weight coefficients from solving the 
secular equation for the sample. Since this is extremely fast, the above method suits for 
simulating very large systems. However, one has to keep in mind that the validity of the 
above relationship relies on imaging the frontier orbital with a bias voltage including only 
the corresponding molecular orbital, and that the STM tip is close enough to the molecule 
that the current flows mostly through the molecular state. 

It is possible to get quite similar dependence for the Tersoff-Hamann approximation for 
simulating STM images. In the TH scheme surfaces of constant local density of sample 
states (LDOS) at the position of the tip apex can be identifled with the STM images. In 
this case the LDOS is dominated by u', so the current will be proportional to {rtlu'). While 



this bears close resemblance to Eq. 2.7 there may be crucial differences especially close to 



the surface, as discussed in Sec. IV B 



III. MODEL AND METHODS 



In this work, simulations have been carried out using our purpose-built method for simu- 
lating STM in a TB atomic orbital basis, and compared to results obtained using standard 
DFT methods (VASP cod^^. In this section we flrst describe the structures of the studied 
systems and how they have been obtained, then explain the methods used to calculate the 
electronic structure of the system - including the molecular orbitals - and flnally detail the 
methods used to calculate the STM images. 

The studied system consists of a Cu(lOO) substrate with a bilayer of NaCl upon which a 
6-gold-pentacene complex is adsorbed. The used geometry for the system has been obtained 
through DFT calculation^. The super cell includes a slab of a four layers of copper each 
layer having 6x9 atoms, and two layers of NaCl each layer consisting of 4x6 Na and CI 
ions. In the studied geometry, seen in Fig. [T| the gold atom is connected to the central 



6 



ring of the pentacene molecule, with the hydrogen tilting upwards out of its way. The 
plane of pentacene lies 3.5 A above the insulating overlayer, and 9 A from the topmost 
copper atoms. The same geometry has been used for calculations of isolated molecules as 
the adsorbed geometry differs only slightly from the relaxed isolated molecule. In the TB 
calculations the STM tip has been simulated using a pyramidal structure consisting of five 
copper atoms, connected to a slab of two layers of Cu(lOO) surface. 

(a)| 




Figure 1: (Color Online) The geometric structure used in the TB simulations of 6-gold-pentacene 
on a NaCl -surface from two directions. Na, CI, Cu, C, H and Au atoms are shown as blue, green, 
brown, grey, white and yellow spheres, respectively. 



Electronic structure calculations for the system, utilizing TB formalism, have been done 
using angle- and distance-dependent Slater-Koster typd^ hopping integrals Vajjm- The hop- 
ping integrals have been modified by including an exponential damping to improve repre- 
sentation of interactions in the insulating overlayer and the molecule. 
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Here a and (3 refer to atomic orbitals on atoms a and b, m indicates the type of bond (a, tt), tq 
designates the distance at which the hopping integrals are changed to behave exponentially, 
A the speed at which the change happens and k the strength of the exponential damping. 
In this study the values for the parameters have been: ro = 3 A, k = 1.0 and A = 2.8 
A^^. A cut-off radius for the hopping integral was used to fasten up the simulations. 

The bias voltage V used in the STM experiments defines the energy range [i^p, -^f + eV^] 
and the states that can contribute to the STM image, and thus restricts the basis set needed 
in the modelling. In all TB calculations, for C and CI valence electrons of s- and p-types have 
been included while for Cu, Au, H and Na only valence electrons of s-type are accounted for. 
Au d-electrons have been neglected since the DFT simulations indicate that the d-orbitals 



of Au do not have a significant contribution near tfie Fermi energy of the system and do 
not therefore affect the STM imaging with reasonably small bias. Similarly for Cu and Na 
only the s-orbitals have contributions near the Fermi energy. The tip apex Cu-atom has also 
been treated as having an s-type orbital only. 

In the Green's function approach, the electronic structure for the Cu atoms in the two 
lowest layers of the adsorbate and for all the atoms in tip is modelled using Haydock's 
recursion schemd^'^ which enables description of a larger background and also speeds up 
the simulations. The two topmost layers of Cu atoms in the adsorbate are modelled similarly 
to the rest of system. 

The TB calculations are carried out with a very small number of parameters — in addition 
to fixed parameters describing the hopping integrals, the on-site energies are required for 
obtaining the Hamiltonian. In general, the on-site energies are not directly fitted to any 
external data but first taken from literatur^^^l a^^j then shifted self-consistently to keep the 
partial charge of each ion neutral. This procedure is used in calculation of the free molecules. 

However, this method does not work for systems with charge transfers between ions, and 
the ionic charges have to be obtained from DFT or experiments. The on-site energies for 
charged Na and CI ions were thus adjusted to reproduce the energy-dependent density of 
states projected to atomic orbitals obtained from DFT calculations. The fitting was done 
only for clean NaCl bilayer on copper. The partial charges for each Na and CI ion were then 
obtained by integrating the partial density of states. 

In calculation of the combined adsorbate-substrate system the on-site energies are again 
shifted self-consistently to maintain the predefined partial charges (zero for all other ions 
except the ions in the polar film, for which the values obtained for the clean film is used). 
Thus there are actually no parameters to be fitted for the combined system, which increases 
the prediction power of the method. 

Previously, the functional forms of the hopping integrals used in the TB method had 
been tested for the clean NaCl bilayer and free organic molecules. The obtained molecular 
orbitals accurately match DFT Kohn-Sham states in shape and qualitatively in energy for 
chloronitrobenzen^^ and pentacen^S!_ Finally, because DFT calculations with and without 
spin-polarization gave similar results for STM simulations, the TB calculations were done 
without taking spin into account^^. 

The calculation of the tunneling current using the TB method is done according to the 
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Todorov-Pendry approach. For comparison, the STM calculations using TB have been 
done using two kinds of functional forms for the hopping integral between the tip and the 
sample. The first is the same modified Slater-Koster typd^ hopping which is used in the 
electronic structure calculations, with exponential damping. The second method is through 
the Wolfsberg-Helmholz approximation using overlap integrals for Slater-type atomic orbitals 
(STO), similarly to extended Hiickel theory (EHT).'^^'^ These hopping integrals have the 
form: 

^a/3(^afe) = ^kSalB {H^a + Hp^) , (3.2) 

where k = 1.75. The overlaps Sa/3 are calculated using the well-known analytical forms 
presented by MuUiken et al!^. The Clementi-Raimondi^ screening constants were used for 
all species of atoms. The H^a and Hp/3 elements are the ionization energies for the atomic 
orbitals in question. 

IV. RESULTS 

A. Molecular Orbital Decompositions 

Molecular orbital calculations can provide a powerful tool in the analysis of STM results 
on insulating overlayers since the conductivity of molecules adsorbed on such surfaces is 
usually heavily influenced by the frontier molecular orbital^. The gold-pentacene complex 
has a singly occupied molecular orbital (SOMO) which dominates the STM current on 
small biases. In Ref. [3] the SOMO state was argued to have significant contributions from 
the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital 
(LUMO) of pentacene, in addition to the gold 6s-orbital. 

In Table |l| can be seen the weight coefficients of free pentacene molecular orbitals, and 
the gold 6s-orbital to the SOMO state of the Au-pentacene complex. Interestingly it can 
be seen that the SOMO state is formed as a linear combination of both the HOMO and 
LUMO orbitals with almost equal weights, with the gold s-orbital also contributing greatly. 
This also verifies the assumptions of Ref. |3| made using geometrical arguments. From the 
square sum of these coefficients (0.91) can be deduced that the other pentacene states that 
lie near the Fermi energy, have only a very small contribution to the SOMO. The states 
above and below in energy to the SOMO (designated here as SOMO-1 and SOMO+1) can 
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be seen to originate from the pentacene LUMO and HOMO states, respectively, with a large 
contribution from the gold 6s orbital and a smaller one from other pentacene states. 





SOMO-1 


SOMO 


SOMO+1 


HOMO 


0.69 


0.61 


0.20 


LUMO 


0.25 


0.66 


0.59 


Au6s 


0.34 


0.31 


0.62 



Table I: The contributions of the HOMO and LUMO of the free pentacene molecule and the 6s 
orbital of the gold atom to the SOMO state which is formed as Au bonds with pentacene. The 



values have been calculated according to Eq. 2.1 in Sec. Uli 



The TB and DFT calculations for SOMO molecular orbital of a free gold-pentacene 
system (as seen in Fig. [2] (a) and (b)) are found to be in good accordance with each other, 
testifying of a successful fitting for the electronic structure in the TB method. The only 
notable difference between the results is that the gold s-orbital is distinctly contracted in 
the TB calculations compared to the DFT results. In the free molecule calculations there 
is a bridge combining two orbital lobes with equal signs, clearly seen in the center of Figs. 
|2ja) and (b). The bridge runs along the long axis of the molecule, seen vertically in the 
pictures. 

In Fig. |2](c) and (d) is presented the molecular orbital calculations for the corresponding 
SOMO* state of the system consisting of the adsorbed molecule, the insulating overlayer 
and the copper substrate. Similarly to the free molecule calculations, the results are in good 
accordance with each other except that the gold s-orbital still being distinctly contracted 
in the TB calculations compared to the DFT. In the adsorption, the bridge seen in the 
free molecule calculations has disappeared entirely, leaving a cross-shaped nodal area. The 
bridge is also absent in experimental STM images of the SOMO* state. 

This absence of the bridge happens due to interaction of the molecule with the chlorine 
atom which lies on the axis connecting gold with its nearest carbon. Simulations done using 
both VASP and TB (with the chlorine interacting solely with Au) using a system consisting 
of the chlorine in question and the Au-pentacene complex exhibit the same behaviour of the 
bond breaking, as the atom is brought closer to the molecule from further away. A similar 
shift in energy is found for the chlorine p-orbital using either DFT or TB with predefined 
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(b) 







Figure 2: (Color Online) Isosurfaces of molecular orbitals of free Au-pentacene complex (SOMO) 
calculated using (a) TB and (b) DFT, and adsorbed Au-pentacene complex (SOMO*) with (c) TB 
and (d) DFT. In the TB molecular orbitals the lobes with different signs have been plotted using 
red and white. The molecular orbitals in the DFT calculations correspond to single Kohn-Sham 
states. 

charge maintenance. In effect, this can be understood as a sign that some interactions 
between adsorbed molecular systems and the insulating overlayer can affect the molecular 
states, and therefore also the STM images which are integrally connected to them. 



B. STM simulations 



The constant current STM images showing the SOMO* orbital simulated using TP ap- 
proach of Eq, 2.2 are shown in Figs. [3|^a) and (b) calculated with the two different hopping 



integrals between the STM tip and the adsorbate. For comparison, in Fig. |3]^c) is shown 
the corresponding simulations carried out using DFT within TH approximation. The STM 
images calculated with DFT and STO-EHT hopping integrals can be seen to be similar 
to each other, while the image calculated with the Slater-Koster hopping integrals clearly 
differs from them — especially near the bridge area. The DFT and STO-EHT images also 
clearly resemble the experimental STM result^ and the molecular orbital calculated for the 
adsorbed molecule (see Fig. |2|. 

Figure |3] (d) presents a plotting of the hopping map to the SOMO* orbital using STO- 
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Figure 3: (Color Online) Simulated STM images of the 6-gold-pentacene complex calculated using 
(a) TB-TP with Slater-Koster type hopping integrals, (b) TB-TP with STO-EHT type hopping 
integrals and (c) DFT-TH. In (d) is shown a hopping map of the adsorbed molecular orbital 
calculated using STO-EHT hopping integrals. In all images the height corresponds to the distance 
from the topmost insulating overlayer. In (d) the color axis is extended in two directions to visualize 
the phase information. The background in the hopping map has been set to zero to clarify the 
image. 

EHT hopping integrals. The presented image corresponds to a constant current STM image. 
To clarify the revealed phase information, the sign of the hopping integral is visualized by 
extending the color axis in two directions. The background, which would go to minus infinity, 
is set to zero to help interpreting the image. The presented hopping map is nearly identical 
to the corresponding STM image, due to the SOMO* state dominating the STM image, and 
it is of the order of 100 times faster to calculate than the TP image of Fig. [sj^b). 

Both TP with STO-EHT and the hopping maps reproduce nicely the arc of the exper- 
imental STM images around the Au atom. Similarly to the molecular orbitals, presented 
in Fig. |2j the arc originating from the Au 6s-orbital is slightly smaller in the TB calcula- 
tions compared to DFT. This is due to either the too weak interaction strength between the 
Au s-orbital and the molecule or too small extension of the the Au Slater orbital used are 
calculated in the electronic structure calculations. Notice that no especial fitting was done 
for the calculation but the parameters for Au are directly taken from literature. The arc is 
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visible in the TB images only close to the surface and thus the TB and DFT calculations 
are done at slightly different distances. 

STM images for TB have been calculated with bias 0.08 V and for DFT with 0.5 V^. 
One should note, that the theoretical values are not directly comparable with experimental 
voltages, as the Coulomb blockade in the experiments shifts the observed energies of the 
states^^. The bias dependence of the STM image is small for both modelling methods with 
reasonably low biases. The corrugations of the images can be seen to be very similar in 
all images. The TB images differ in size from the DFT images, due to the lower imaging 
distance. 

Even though there is little bias dependence for the STM image of the system, the image 
may be quite dependent on the chosen constant current surface, i.e. the initial height from 
the adsorbate. To visualize this we have calculated the STM current in the plane where the 
C-Au bond lies. The result is shown in Fig. |4]as a logarithmic mapping of the STM current 
as a function of height. There are several interesting properties to the STM image, which 
can be seen from this figure. The most striking property is the fact that the STM current is 
clearly not unambiguous; there exist folded constant current surfaces — a predefined current 
does not give a unique height. This might have implications in experiments, such as that 
the scanning direction affects the outcome of the measurement. 
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Figure 4: (Color Online) Logarithmic plotting of the STM current as a function of height in the 
plane of the Au-Pentacene bond, calculated using TB with STO-EHT hopping integrals. Positions 
of carbons in the plane marked with crosses, hydrogens with circles and gold with a square. Lines 
correspond to constant current surfaces. An ambiguous surface is highlighted. 

One should note, that according to the geometry obtained in the DFT calculations, in 
this system the highest atom lies 4.5 A above the insulating overlayer, and approximately 10 
A from the atoms in the conducting surface. Because of this, the distance between the STM 
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tip and the molecule cannot be very large while still being able to generate a measurable 
tunneling current. As can be seen in Fig. |4| an ambiguity can be found also quite far away 
from the molecule itself. Thus, it is possible that the experiments are done within a range 
where the above multiplicity can occur. However, being close to the molecule implies that 
the actual effect on experiments will depend on the real structure of the STM tip. 

In addition to this, the contrast and extension of the STM image clearly depend on the 
followed constant current surface. As is well-known, the STM image is seen to widen further 
away from the molecule, and the contrast of the nodal region in the middle vanishes at higher 
distances. However, due to complicated shape of the SOMO orbital, the topographical 
STM image may now be also qualitatively different with various values of the tunneling 
current. Thus, this may give an answer to the experimentally observed phenomena of the 
arc appearing in some experimental STM images while being absent in others. From the 
logarithmic mapping, as well as the calculated STM images with different initial heights, 
can be interpreted that the arc should be seen at low scanning heights but disappear when 
going further away from the surface. 

The STM images calculated using the two hopping integrals (Fig. [3](a) and (b)) empha- 
size the importance of the proper choice of the hopping integral in TB calculations. The 
differences between the two approximations can be further visualized with the hopping maps 
whose cross sections along the plane of the Au-C bond are shown in Fig. [s] (a) and (b) for 
the used Slater-Koster and STO-EHT type hopping integrals, respectively. The hopping 
map results differ drastically. Where in the STO-EHT a slight negative protrusion rises 
slightly between the two positive lobes, in the second image it is a large negative lobe that 
encases the positive lobe by connecting to the negative lobe centered around the gold atom. 
This is the reason for the clear differences between the STM images for the two methods in 
Fig. [3j This is partly due to the incorrect spatial dependence of the Slater-Koster hopping 
integrals neglecting the size (principal quantum number) of the orbitals. This is especially 
important in this kind of a system in which a heavy atom (gold and its 6s-orbital) plays 
a major role. On a brighter note, for this kind of systems the STO-EHT hopping integral 
gives a good accordance with both DFT and experimental results. 

For comparison, a similar cross-section of the SOMO* state is shown in Fig. ^c). Even 
though similarity with STO-EHT hopping map further away from the surface, there are 
some discrepancies close to the molecule. This highlights the difference in Tersoff-Hamann 
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Figure 5: (Color Online) Images of the spatial dependence of the (a) Slater-Koster type hopping 
integrals to the SOMO* orbital, (b) STO-EHT type hopping integrals to the SOMO* orbital and 
(c) the local density of the SOMO* orbital (calculated using TB) as a function of height in the 
plane of the Au-pentacene bond. Locations of carbons in the plane marked with crosses, hydrogens 
with circles and gold with a square. The values of the hopping integrals and the molecular orbital 
have been truncated to help interpretation of the image. 

approach compared to Todorov-Pendry approximation. 

By modifying the range at which the exponential tail is applied and by changing the 
relative differences of s-p and s-s interactions the Slater-Koster type hopping integrals can 
be brought into line with the DFT and STO-EHT results, producing similar STM images and 
hopping maps. However, this greatly affects the ability of the method to predict measured 
STM images, especially compared to the STO-EHT, for which no specific fitting was carried 
out. This emphasizes the importance of a proper modelling of the tip-adsorbate interaction. 
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V. SUMMARY 



The objective of this work has been to develop TB based methods for simulating STM 
images of complicated molecules on insulating overlayers. We used pentacene bonded with 
gold on a sodium chloride thin film with copper substrate as a test system. Results obtained 
using our purpose-built method for simulating STM images are found to be sensitive to the 
choice of hopping integral between the STM tip and the adsorbate. 

With STO-EHT -type hopping integrals the STM results very closely match both exper- 
iments and DFT simulations. The only notable difference found is regarding the size and 
visibility of the arc formed by the Au 6s-orbital. Instead, the Slater-Koster -type hopping 
integrals give clearly differing results to both other simulations and experiments, due to their 
disregard for the size of the orbitals. 

A logarithmic plotting of the STM current with the STO-EHT hopping integrals predicts 
that the complex structure of the molecular orbitals may lead to ambiguity of the STM 
image. This kind of effect may have implications in experiments, such as the STM image 
changing with respect to the scanning direction. Additionally, the tip may drift too close 
to the molecule in constant current imaging, allowing the molecule to jump to the tip. The 
actual effects of the ambiguity will largely depend on the actual shape of the tip, which would 
merit a further study. These results highlight the importance of combining experiments with 
simulations in order to properly interpret the obtained results. 

Hopping maps, for which the theory has been presented in this work, provide a powerful 
tool, not only in their abihty to study the phase of the STM current, but also in the fact 
that they can be used to very easily simulate the STM results of different conformations 
even for large systems. As can be seen from the results presented here, for systems where a 
single molecular orbital dominates the STM image, the hopping map is very similar to the 
STM image. The phase of the STM current can be utilized to analyze the formation of the 
STM current from different parts. 

The methods developed in this work will be used to study large molecules suitable for 
molecular electronics on insulating overlayers. The ability to simulate large systems with 
comparably small computational resources, especially when utilizing the hopping maps, 
makes it more feasible to simulate relevant single-molecular systems. 
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Appendix A: DIFFERENT FORMULATIONS OF TODOROV-PENDRY AP- 
PROACH 



The current running through a system is obtained in Todorov-Pendrj^ approach from 



2716 

IT 



t-Ep+eV 

/ ^ApIAE - eV)Tl,^,pl,^{E)T,;\dE. (Al) 

JEf 



In tight-binding basis, the density matrix p is easily calculated with the aid of Green's 
function, 

TT ^ — ' ■' 

s 

where i, f denote initial (for example, in the tip) and final state (in the sample). 



With the aid of Eqs. Al and A2 it is easy to transform the integrand into form 



J = ^ np'r'rT^pL'T^'r'] = ^ Tr [p°,,T4G°7G5+ T.^'] • (A3) 
Various other formulations to calculate tunneling current can be recovered from the equa- 



tion above. If we apply K2] once again, we get the following form: 



J = ^Tr[Gt^Tt^G°7Gj^,T.,.,G^;:] 
Assuming that tip and sample are disconnected in the beginning and utilizing Dyson's 



equation: 



since there are no tip-sample matrix elements in the case of uncoupled system. 



Thus it is straightforward to convert Eq. A3 to the form 



J = ^ np'r'rVlG-jG%y^.^,], (A4) 

which is essentially the formalism used in this work. 

There are some extensions and alternative formulations that can be shown here. First, the 
convergence parameter t] can be generalized to an energy dependent broadening parameter 
r/2 related to self-energy. Hence we can express this in a Landauer-Biittiker form 

J = I np\V}.G-f,Tf,fG%yM (A5) 
as formulated by Meir and WingreenP^l. 
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In addition, applying Dyson's equation leads to the equation 

which is the result given by Ness and Fishei'^. 

One of the equivalent formulations is given by McKinnon and Choj^, which has been 
written in a more elegant form by Mingo et al?^ . A further application of Dyson's equation 



to Eq. |A3| tells us 

G-j = G^^iV^i jiG^j, J 

— r^o- I T/t r^o- 
^i'i ~ ^i'i ' '^if'^ f"i"'^i"i 

which means 

Gii = G% (/ - V^^,j,G'j,j,Vj„,-„G%^,) ^ 

and 

r^- _ r<0- ( T _ T/t r^O- T/t 1^0- "i-lT/t /^O- _ i^O~rp-\ ^0- 
Following Mingo et al^ we utilize a shorthand notation 

T't — iT _ T/t r'^^ T/t r''^^ "i^^T/t — n~ T/t 

Which leads to the form: 

j = ^ Tr[p°,,D;,r4G°7G;+ = ^ Tr[pl,^D;X.PL'DX'V^'r']- (A6) 

Now it mostly depends on the computational method, which one is the most suitable 
formalism to calculate the tunneling current. 
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